#include "stdinc.h"
#include "newhash.h"
#include "extfunc.h"
#include "extvab.h"

static char src_rc_seq[1024];
extern long long single_count;
extern long long single_map;
void readseq1by1(char *src_seq, char *src_name, int *len_seq, FILE *fp,long long num_seq)
{
	int	i,k, n,strL;
	char	c;
	char	str[5000];

	n = 0;
	k = num_seq;
	while(fgets(str, 4950, fp))	{
		if(str[0] == '#')	continue;
		if(str[0] == '>')	{
			/*
			if(k >= 0)	{  // if this isn't the first '>' in the file
				*len_seq = n;
			}
			*/
			*len_seq = n;
			n = 0;
			sscanf(&str[1],"%s",src_name);
			return;
		} else {
			strL = strlen(str);
			if(strL+n>maxReadLen)
				strL = maxReadLen - n;
			for(i = 0; i < strL; i ++)	{
				if(str[i] >= 'a' && str[i] <= 'z') {
					c = base2int(str[i]-'a'+'A');
					src_seq[n ++] = c;
				} else if(str[i] >= 'A' && str[i] <= 'Z') {
					c = base2int(str[i]);
					src_seq[n ++] = c;
				          // after pre-process all the symbles would be a,g,c,t,n in lower or upper case.
				} else if(str[i]=='.') {
					c = base2int('A');
					src_seq[n ++] = c;
				}          // after pre-process all the symbles would be a,g,c,t,n in lower or upper case.
			}
			//printf("%d: %d\n",k,n);
		}
	}

	if(k >= 0){
		*len_seq = n;
		return;
	}
	*len_seq = 0;
}


void read_one_sequence(FILE *fp, long long *T, char **X)

{

	char *fasta,*src_name; //point to fasta array
	int num_seq,len,name_len,min_len;

	num_seq = readseqpar(&len,&min_len,&name_len,fp);
	if(num_seq<1){
		printf("no fasta sequence in file\n");
		*T = 0;
		return;
	}
	fasta = (char *)ckalloc(len*sizeof(char));
	src_name = (char *)ckalloc((name_len+1)*sizeof(char));
	rewind(fp);

	readseq1by1(fasta,src_name,&len,fp,-1);
	readseq1by1(fasta,src_name,&len,fp,0);
	
	*X = fasta;
	*T = len;
	free((void *)src_name);
}

long long multiFileParse(int *max_leg, int *min_leg,int *max_name_leg, FILE *fp)
{
	
	char str[5000];
	FILE *freads;
	int slen;
	long long counter = 0;
	*max_name_leg = *max_leg = 1;
	*min_leg = 1000;
	while(fgets(str,4950,fp)){
		slen = strlen(str);
		str[slen-1] = str[slen];
		freads = ckopen(str,"r");
		counter += readseqpar(max_leg,min_leg,max_name_leg,freads);
		fclose(freads);
	}
	return counter;
}

long long readseqpar(int *max_leg, int *min_leg,int *max_name_leg, FILE *fp)
{
	int	l, n;
	long long k;
	char	str[5000], src_name[5000];


	n = 0;
	k = -1;
	while(fgets(str, 4950, fp))	{
		if(str[0] == '>')	{
			if(k >= 0)	{
				if(n > *max_leg)	
					*max_leg = n;
				if(n < *min_leg)	
					*min_leg = n;
			
			}
			n = 0;
			k ++;
			sscanf(&str[1], "%s", src_name);
			if((l = strlen(src_name)) > *max_name_leg)		
				*max_name_leg = l;
		} else {
			n += strlen(str)-1;
		}
	}
	if(n > *max_leg)	
		*max_leg = n;
	
	if(n < *min_leg)	
		*min_leg = n;

	k ++;
	return(k);
}

void read1seqfq(char *src_seq, char *src_name, int *len_seq, FILE *fp)
{
	int	i,n,strL;
	char	c;
	char	str[5000];
	boolean flag=0;

	while(fgets(str, 4950, fp))	{
		if(str[0]=='@'){
			flag = 1;
			sscanf(&str[1],"%s",src_name);
			break;
		}
	}

	if(!flag){   //last time reading fq file get this
		*len_seq = 0;
		return;
	}

	n = 0;
	while(fgets(str, 4950, fp)){
		if(str[0] == '+')	{
			fgets(str,4950,fp);   // pass quality value line
			*len_seq = n;
			return;
		} else {
			strL = strlen(str);
			if(strL+n>maxReadLen)
				strL = maxReadLen - n;
			for(i = 0; i < strL; i ++)	{
				if(str[i] >= 'a' && str[i] <= 'z') {
					c = base2int(str[i]-'a'+'A');
					src_seq[n ++] = c;
				} else if(str[i] >= 'A' && str[i] <= 'Z') {
					c = base2int(str[i]);
					src_seq[n ++] = c;
				          // after pre-process all the symbles would be a,g,c,t,n in lower or upper case.
				} else if(str[i]=='.') {
					c = base2int('A');
					src_seq[n ++] = c;
				}          // after pre-process all the symbles would be a,g,c,t,n in lower or upper case.
			}
			//printf("%d: %d\n",k,n);
		}
	}

	*len_seq = n;
	return;
}

// find the next file to open in libs
static int nextValidIndex(int libNo,boolean pair,unsigned char asm_ctg)
{
	int i=libNo;
	
	while(i<num_libs){
		if(asm_ctg==1&&(lib_array[i].asm_flag!=1&&lib_array[i].asm_flag!=3)){
			i++;
			continue;
		}else if(asm_ctg==0&&(lib_array[i].asm_flag!=2&&lib_array[i].asm_flag!=3)){
			i++;
			continue;
		}else if(asm_ctg>1&&lib_array[i].asm_flag!=asm_ctg){  // reads for other purpose 
			i++;
			continue;
		}
		if(lib_array[i].curr_type==1&&
			lib_array[i].curr_index<lib_array[i].num_a1_file)
			return i;
		if(lib_array[i].curr_type==2&&
			lib_array[i].curr_index<lib_array[i].num_q1_file)
			return i;
		if(lib_array[i].curr_type==3&&
			lib_array[i].curr_index<lib_array[i].num_p_file)
			return i;
		if(pair){
			if(lib_array[i].curr_type<3){
				lib_array[i].curr_type++;
				lib_array[i].curr_index = 0;
			}else
				i++;
			continue;
		}
		if(lib_array[i].curr_type==4&&
			lib_array[i].curr_index<lib_array[i].num_s_a_file)
			return i;
		if(lib_array[i].curr_type==5&&
			lib_array[i].curr_index<lib_array[i].num_s_q_file)
			return i;
		if(lib_array[i].curr_type<5){
			lib_array[i].curr_type++;
			lib_array[i].curr_index = 0;
		}else
			i++;
	}//for each lib
		
	return i;
}

static FILE *openFile4read(char *fname)
{
	FILE *fp;
	if(strlen(fname)>3&&strcmp(fname+strlen(fname)-3,".gz")==0){
		char *cmd = (char *)ckalloc((strlen(fname)+20)*sizeof(char));
		sprintf(cmd,"gzip -dc %s",fname);
		fp = popen(cmd,"r");
		free(cmd);
		return fp;
	}else{
		return ckopen(fname,"r");	
	}
	
}

void openFileInLib(int libNo)
{
	int i = libNo;
	if(lib_array[i].curr_type==1){
		printf("[%s]opened file:\n %s\n",
					__FUNCTION__,lib_array[i].a1_fname[lib_array[i].curr_index]);
		printf("[%s]opened file:\n %s\n",
					__FUNCTION__,lib_array[i].a2_fname[lib_array[i].curr_index]);
		lib_array[i].fp1 = openFile4read(lib_array[i].a1_fname[lib_array[i].curr_index]);	
		lib_array[i].fp2 = openFile4read(lib_array[i].a2_fname[lib_array[i].curr_index]);	
		lib_array[i].curr_index++;
		lib_array[i].paired = 1;
	}else if(lib_array[i].curr_type==2){
		printf("[%s]opened file:\n %s\n",
					__FUNCTION__,lib_array[i].q1_fname[lib_array[i].curr_index]);
		printf("[%s]opened file:\n %s\n",
					__FUNCTION__,lib_array[i].q2_fname[lib_array[i].curr_index]);
		lib_array[i].fp1 = openFile4read(lib_array[i].q1_fname[lib_array[i].curr_index]);	
		lib_array[i].fp2 = openFile4read(lib_array[i].q2_fname[lib_array[i].curr_index]);	
		lib_array[i].curr_index++;
		lib_array[i].paired = 1;
	}else if(lib_array[i].curr_type==3){
		printf("[%s]opened file:\n %s\n",
                       __FUNCTION__,lib_array[i].p_fname[lib_array[i].curr_index]);
		lib_array[i].fp1 = openFile4read(lib_array[i].p_fname[lib_array[i].curr_index]);	
		lib_array[i].curr_index++;
		lib_array[i].paired = 0;
	}else if(lib_array[i].curr_type==4){
		printf("[%s]opened file:\n %s\n",
					__FUNCTION__,lib_array[i].s_a_fname[lib_array[i].curr_index]);
		lib_array[i].fp1 = openFile4read(lib_array[i].s_a_fname[lib_array[i].curr_index]);	
		lib_array[i].curr_index++;
		lib_array[i].paired = 0;
	}else if(lib_array[i].curr_type==5){
		printf("[%s]opened file:\n %s\n",
					__FUNCTION__,lib_array[i].s_q_fname[lib_array[i].curr_index]);
		lib_array[i].fp1 = openFile4read(lib_array[i].s_q_fname[lib_array[i].curr_index]);	
		lib_array[i].curr_index++;
		lib_array[i].paired = 0;
	}
	
}

static void reverse2k(char *src_seq,int len_seq)
{
	if(!len_seq)
		return;
	
	int i;
	reverseComplementSeq(src_seq,len_seq,src_rc_seq);

	for(i=0;i<len_seq;i++)
		src_seq[i] = src_rc_seq[i];
}

static void closeFp1InLab(int libNo)
{
	int ftype = lib_array[libNo].curr_type;
	int index = lib_array[libNo].curr_index-1;
	char *fname;
	if(ftype==1)
		fname = lib_array[libNo].a1_fname[index];
	else if(ftype==2)
		fname = lib_array[libNo].q1_fname[index];
	else if(ftype==3)
		fname = lib_array[libNo].p_fname[index];
	else if(ftype==4)
		fname = lib_array[libNo].s_a_fname[index];
	else if(ftype==5)
		fname = lib_array[libNo].s_q_fname[index];
	else
		return;
	if(strlen(fname)>3&&strcmp(fname+strlen(fname)-3,".gz")==0)
		pclose(lib_array[libNo].fp1);
	else
		fclose(lib_array[libNo].fp1);
}

static void closeFp2InLab(int libNo)
{
	int ftype = lib_array[libNo].curr_type;
	int index = lib_array[libNo].curr_index-1;
	char *fname;
	if(ftype==1)
		fname = lib_array[libNo].a2_fname[index];
	else if(ftype==2)
		fname = lib_array[libNo].q2_fname[index];
	else
		return;
	if(strlen(fname)>3&&strcmp(fname+strlen(fname)-3,".gz")==0)
		pclose(lib_array[libNo].fp2);
	else
		fclose(lib_array[libNo].fp2);
}

boolean read1seqInLib(char *src_seq, char *src_name, int *len_seq, int *libNo,boolean pair,unsigned char asm_ctg)
{
	int i = *libNo;
	int prevLib = i;

	if(!lib_array[i].fp1  // file1 does not exist 
		||(lib_array[i].curr_type!=1&&feof(lib_array[i].fp1))   // file1 reaches end and not type1
			||(lib_array[i].curr_type==1&&feof(lib_array[i].fp1)&&feof(lib_array[i].fp2))){//f1&f2 reaches end
		if(lib_array[i].fp1&&feof(lib_array[i].fp1)){
				 closeFp1InLab(i);
				 //printf("[%s]%d reads in current file , (%.1f) map-rate .\n",__FUNCTION__,single_count,single_map/single_count);
				 single_count=single_map=0;
		}
		if(lib_array[i].fp2&&feof(lib_array[i].fp2)){
				closeFp2InLab(i);
				//printf("[%s]%d reads in current file , (%.1f) map-rate .\n",__FUNCTION__,single_count,single_map/single_count);
				single_count=single_map=0;
		}
		
		*libNo = nextValidIndex(i,pair,asm_ctg);
		i = *libNo;
		if(lib_array[i].rd_len_cutoff>0)
			maxReadLen = lib_array[i].rd_len_cutoff<maxReadLen4all ?
								lib_array[i].rd_len_cutoff:maxReadLen4all;
		else
			maxReadLen = maxReadLen4all;

		//record insert size info
		//printf("from lib %d to %d, read %lld to %ld\n",prevLib,i,readNumBack,n_solexa);
		if(pair&&i!=prevLib){
			if(readNumBack<n_solexa){
				pes[gradsCounter].PE_bound = n_solexa;
				pes[gradsCounter].rank = lib_array[prevLib].rank;
				pes[gradsCounter].pair_num_cut = lib_array[prevLib].pair_num_cut;
				pes[gradsCounter++].insertS = lib_array[prevLib].avg_ins;
				readNumBack = n_solexa;
			}		
		}
		if(i>=num_libs)
			return 0;
		openFileInLib(i);

		if(lib_array[i].curr_type==1){
			readseq1by1(src_seq, src_name,len_seq, lib_array[i].fp1,-1);
			readseq1by1(src_seq, src_name,len_seq, lib_array[i].fp2,-1);
		}else if(lib_array[i].curr_type==3||lib_array[i].curr_type==4)
			readseq1by1(src_seq, src_name,len_seq, lib_array[i].fp1,-1);

	}
	if(lib_array[i].curr_type==1){
		if(lib_array[i].paired==1){
			readseq1by1(src_seq, src_name,len_seq, lib_array[i].fp1,1);
			if(lib_array[i].reverse)
				reverse2k(src_seq,*len_seq);
			lib_array[i].paired = 2;
			if(*len_seq>0||!feof(lib_array[i].fp1)){
				n_solexa++;
				return 1;
			}
			else
				return read1seqInLib(src_seq,src_name,len_seq,libNo,pair,asm_ctg);
		}else{
			readseq1by1(src_seq, src_name,len_seq, lib_array[i].fp2,1);

			if(lib_array[i].reverse)
				reverse2k(src_seq,*len_seq);
			lib_array[i].paired = 1;
			n_solexa++;
			return 1; //can't fail to read a read2
		}
	}
	if(lib_array[i].curr_type==2){
		if(lib_array[i].paired==1){
			read1seqfq(src_seq, src_name,len_seq, lib_array[i].fp1);
			/*
			if(*len_seq>0){
				for(j=0;j<*len_seq;j++)
					printf("%c",int2base(src_seq[j]));	
				printf("\n");
			}
			*/
			if(lib_array[i].reverse)
				reverse2k(src_seq,*len_seq);
			lib_array[i].paired = 2;
			if(*len_seq>0||!feof(lib_array[i].fp1)){
				n_solexa++;
				return 1;
			}else
				return read1seqInLib(src_seq,src_name,len_seq,libNo,pair,asm_ctg);
		}else{
			read1seqfq(src_seq, src_name,len_seq, lib_array[i].fp2);
			if(lib_array[i].reverse)
				reverse2k(src_seq,*len_seq);
			lib_array[i].paired = 1;
			n_solexa++;
			return 1; //can't fail to read a read2
		}
	}
	if(lib_array[i].curr_type==5)
		read1seqfq(src_seq, src_name,len_seq, lib_array[i].fp1);
	else{
		readseq1by1(src_seq, src_name,len_seq, lib_array[i].fp1,1);
	}
	/*
	int t;
	for(t=0;t<*len_seq;t++)
		printf("%d",src_seq[t]);
	printf("\n");
	*/
	if(lib_array[i].reverse)
		reverse2k(src_seq,*len_seq);
	if(*len_seq>0||!feof(lib_array[i].fp1)){
		n_solexa++;
		return 1;
	}else
		return read1seqInLib(src_seq,src_name,len_seq,libNo,pair,asm_ctg);
}
